{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"http://landlab.github.io\"><img style=\"float: left\" src=\"../../landlab_header.png\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using NHDPlus HR Datasets With the Landlab NetworkSedimentTransporter Component \n",
    "\n",
    "\n",
    "This tutorial illustrates how to use USGS NHDPlus HR datasets with the NetworkSedimentTransporter after they have been preprocessed in QGIS. We will use the \"8_Final_Links.shp\" and \"8_Final_Nodes.shp\" files that were created according to the instructions that accompany this notebook.\n",
    "\n",
    "\n",
    "**Part 1**: Here, we demonstrate how to import the files into Landlab, convert the shapefile attributes, and remove flat areas from the topograpghy. \n",
    "\n",
    "**Parts 2 - 5**: These cells are copied nearly verbatim from the Landlab notebook _\"Using the Landlab NetworkSedimentTransporter component starting with a shapefile river network\"_. The only differences are the values of  _timesteps_ (cell 16), _dt_ (cell 16), and _originating_link_ (cell 20). This is to demonstrate how the code in Part 1 allows the files to behave exatly the same as the MethowSubBasin shapefiles used in other notebooks.\n",
    "\n",
    "**NOTE:** This notebook should be run in the _landlab_dev_ environment. A future version of this notebook may be available on the Landlab website; if so, use that version."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the necessary libraries, plus a bit of magic so that we can plot within this notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import os\n",
    "import pathlib\n",
    "import functools\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from numpy import nan\n",
    "from landlab.components import FlowDirectorSteepest, NetworkSedimentTransporter\n",
    "from landlab.data_record import DataRecord\n",
    "from landlab.grid.network import NetworkModelGrid\n",
    "from landlab.plot import graph\n",
    "from landlab.io import read_shapefile\n",
    "from landlab import ExampleData\n",
    "\n",
    "from landlab.plot import plot_network_and_parcels\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Load the Preprocessed NHDPlus HR Files\n",
    "\n",
    "First, we need to load the two shapefiles (and their .dbf files) you created using QGIS (these are the only lines of code you need to manually change in this notebook)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "link_shp = # path to the \"8_Final_Links.shp\" file\n",
    "link_dbf = # path to the \"8_Final_Links.dbf\" file\n",
    "node_shp = # path to the \"8_Final_Nodes.shp\" file\n",
    "node_dbf = # path to the \"8_Final_Nodes.dbf\" file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we initiate the grid, identify the pour point node, convert the attribute formats so Landlab can read them, and convert some units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = read_shapefile(\n",
    "    file = link_shp,\n",
    "    points_shapefile = node_shp,\n",
    "    dbf = link_dbf,\n",
    "    points_dbf = node_dbf,\n",
    "    link_field_conversion={\"TotDASqKm\" : \"drainage_area\", \n",
    "                           \"Slope\" : \"channel_slope\", \n",
    "                           \"LengthKM\" : \"reach_length\", \n",
    "                           \"MaxElevSmo\" : \"topographic__elevation\",\n",
    "                          },\n",
    "    node_field_conversion={\"TotDASqKm\" : \"drainage_area\", \n",
    "                           \"Slope\" : \"channel_slope\", \n",
    "                           \"LengthKM\" : \"reach_length\", \n",
    "                           \"MaxElevSmo\" : \"topographic__elevation\",\n",
    "                          },\n",
    "    threshold=0.01,\n",
    "    )\n",
    "\n",
    "# IDENTIFY POUR POINT NODE ID\n",
    "watershed_pour_point = np.squeeze(np.where(np.isnan(grid.at_node[\"drainage_area\"].astype(np.float))))\n",
    "\n",
    "# GIVE POUR POINT NODE THE SAME ATTRIBUTES (PLUS / MINUS A SMALL AMOUNT) AS THE NODE IMMEDIATELY UPSTREAM ON THE MAIN STEM\n",
    "grid.at_node[\"drainage_area\"][watershed_pour_point] = (np.max(grid.at_node[\"drainage_area\"][grid.adjacent_nodes_at_node[watershed_pour_point]])) + ( (np.max(grid.at_node[\"drainage_area\"][grid.adjacent_nodes_at_node[watershed_pour_point]])) * 0.001 )    # POUR POINT NODE DRAINAGE AREA IS 0.1% LARGER THAN ITS NEAREST MAINSTEM NEIGHBOR NODE\n",
    "grid.at_node[\"topographic__elevation\"][watershed_pour_point] = np.min(grid.at_node[\"topographic__elevation\"][grid.adjacent_nodes_at_node[watershed_pour_point]]) - np.min(grid.at_link[\"channel_slope\"].astype(np.float))                                   # POUR POINT NODE ELEVATION IS ([NEAREST MAINSTEM NEIGHBOR NODE ELEVATION] - [MINIMUM SLOPE IN NHDPLUS DATASET])\n",
    "\n",
    "# CONVERT NODE ATTRIBUTES FROM OBJECTS TO FLOATS\n",
    "grid.at_node[\"drainage_area\"] = grid.at_node[\"drainage_area\"].astype(np.float)\n",
    "grid.at_node[\"channel_slope\"] = grid.at_node[\"channel_slope\"].astype(np.float)\n",
    "grid.at_node[\"reach_length\"] = grid.at_node[\"reach_length\"].astype(np.float)\n",
    "grid.at_node[\"topographic__elevation\"] = grid.at_node[\"topographic__elevation\"].astype(np.float)\n",
    "\n",
    "# CONVERT UNITS\n",
    "grid.at_link[\"drainage_area\"] = grid.at_link[\"drainage_area\"] * 1000000                   # CONVERT km^2 TO m^2\n",
    "grid.at_node[\"drainage_area\"] = grid.at_node[\"drainage_area\"] * 1000000                   # CONVERT km^2 TO m^2\n",
    "grid.at_link[\"reach_length\"] = grid.at_link[\"reach_length\"] * 1000                        # CONVERT km TO m\n",
    "grid.at_node[\"reach_length\"] = grid.at_node[\"reach_length\"] * 1000                        # CONVERT km TO m\n",
    "grid.at_link[\"topographic__elevation\"] = grid.at_link[\"topographic__elevation\"] / 100     # CONVERT cm to m\n",
    "grid.at_node[\"topographic__elevation\"] = grid.at_node[\"topographic__elevation\"] / 100     # CONVERT cm to m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we have to deal with cetain areas of the grid called \"flats\". These occur where two or more adjacent nodes have the same topographic elevation value, creating an effective slope of zero between them. The NHDPlus HR datasets assign these areas a slope value of 1E-5 so that there are not any zeros:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.nanmin(grid.at_node[\"channel_slope\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the smallest slope value is 1E-5 (we use the numpy.nanmin function because the value of the pour point node is \"nan\"). However, the Landlab flow directors use topographic elevation to independantly calculate slope. This can result in zero values that confuse the flow directors. \n",
    "\n",
    "We can solve this by performing a \"carve\". We will loop through each node on the grid and calculate the downstream path from that location. Wherever we encounter a flat, we will adjust the topographic elevation by subtracting a height equal to the length of that reach multiplyed by the minimum slope assigned by the NHDPlus HR datatables (1E-5)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CREATE ARRAY THAT WILL STORE TERMINAL FLAT AREAS FOR EACH NODE\n",
    "downstream_nodes = np.zeros(grid.number_of_nodes)\n",
    "\n",
    "# FOR EACH NODE, IDENTIFY TERMINAL FLAT NODE\n",
    "for i in grid.nodes:\n",
    "    adjacent_nodes = grid.adjacent_nodes_at_node[i]\n",
    "    adjacent_nodes = np.append(adjacent_nodes, i)\n",
    "    adjacent_nodes = np.delete(adjacent_nodes, np.where(adjacent_nodes == -1))\n",
    "    adjacent_nodes_DA = (grid.at_node[\"drainage_area\"][adjacent_nodes])\n",
    "    downstream_id = np.where(adjacent_nodes_DA == np.max(adjacent_nodes_DA))\n",
    "    downstream_id = functools.reduce(lambda sub, elem: sub * 10 + elem, downstream_id)\n",
    "    if np.size(downstream_id) > 1:                 \n",
    "        downstream_nodes[i] = nan\n",
    "    else:                                      \n",
    "        downstream_nodes[i] = adjacent_nodes[downstream_id]\n",
    "downstream_nodes[watershed_pour_point] = nan\n",
    "downstream_nodes = downstream_nodes.astype(np.int)                              \n",
    "\n",
    "# CREATE ARRAY TO STORE NEW ELEVATIONS\n",
    "carve = np.ones(grid.number_of_nodes) * grid.at_node[\"topographic__elevation\"]\n",
    "\n",
    "# CALCULATE CARVE\n",
    "for i in grid.nodes:\n",
    "    current_node = i\n",
    "    q = 1\n",
    "    while q == 1:\n",
    "        downstream_node = downstream_nodes[current_node]\n",
    "        if downstream_node >= 0:                                        \n",
    "            current_node_elev = carve[current_node]\n",
    "            downstream_node_elev = carve[downstream_node]\n",
    "            if downstream_node_elev < current_node_elev:\n",
    "                current_node = downstream_node\n",
    "            if downstream_node_elev >= current_node_elev:\n",
    "                carve[downstream_node] = carve[current_node] - (grid.at_node[\"channel_slope\"][current_node] * grid.at_node[\"reach_length\"][current_node])    \n",
    "            if current_node == watershed_pour_point:\n",
    "                q = 0\n",
    "        if downstream_node < 0:                                     \n",
    "            q = 0\n",
    "\n",
    "# RECORD AMOUNT CARVED\n",
    "carve_delta = grid.at_node[\"topographic__elevation\"] - carve\n",
    "\n",
    "# PERFORM CARVE\n",
    "grid.at_node[\"topographic__elevation\"] = carve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable \"carve_delta\" records how much height was subtrated from each node. We can use it to view some statistics on how much we altered the topography:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.bar(grid.nodes, \n",
    "         carve_delta)\n",
    "plt.xlabel('Node ID')\n",
    "plt.ylabel('Vertical distance carved (m)')\n",
    "\n",
    "print('Total number of nodes carved = ', np.count_nonzero(carve_delta))\n",
    "print('Percent of nodes carved = ', (np.count_nonzero(carve_delta) / grid.number_of_nodes) * 100, '%')\n",
    "print('Maximum distance carved = ', np.max(carve_delta), 'm')\n",
    "print('Mean distance (+/- 1-sigma) carved = ', np.mean(carve_delta[np.nonzero(carve_delta)]), '+/-', np.std(carve_delta[np.nonzero(carve_delta)]),  'm')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the adjustments we made to the topography were very small (no more than a few millimeters at most). This means it is very unlikely that these changes will significantly affect the model behavior. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Review the Grid and Add Parameters \n",
    "\n",
    "Alright, let's see what fields we read in with this shapefile:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.at_link.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.at_node.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great! Looks like we have length (reach length), upstream drainage area (drainage area), x and y verticies of each link/reach (x and y of polyline), and bed elevation (topographic elevation).\n",
    "\n",
    "Note that \"reach_length\" is defined by the user, rather than calculated as the minimum distance between nodes. This accounts for channel sinuosity. In this case, \"reach_length\" could be equivalently calculated as the cumulative distance between verticies defined by x and y of polyline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graph.plot_graph(grid, at=\"node,link\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.number_of_links"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.number_of_nodes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our network consists of 317 links between 318 nodes. In the plot above, X and Y represent the plan-view coordinates of the node locations. \n",
    "\n",
    "Next, we need to populate the grid with the relevant topographic and hydrologic information: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.at_node[\"bedrock__elevation\"] = grid.at_node[\"topographic__elevation\"].copy()\n",
    "\n",
    "grid.at_link[\"channel_width\"] = 1 * np.ones(grid.number_of_links) # m\n",
    "\n",
    "grid.at_link[\"flow_depth\"] = 0.5 * np.ones(grid.number_of_links) # m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We must distinguish between topographic elevation (the top surface of the bed sediment) and bedrock elevation (the surface of the river in the absence of modeled sediment)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Create sediment 'parcels' in a DataRecord\n",
    "\n",
    "We represent sediment in the network as discrete parcels (or packages) of grains of uniform size and characteristics. Each parcel is tracked through the network grid according to sediment transport and stratigraphic constraints. \n",
    "\n",
    "Parcels are tracked using the Landlab [DataRecord](../data_record/DataRecord_tutorial.ipynb).\n",
    "\n",
    "First, let's create arrays with all of the essential sediment parcel variables: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# element_id is the link on which the parcel begins. \n",
    "element_id = np.repeat(np.arange(grid.number_of_links), 50)\n",
    "element_id = np.expand_dims(element_id, axis=1)\n",
    "\n",
    "volume = 1*np.ones(np.shape(element_id))  # (m3)\n",
    "active_layer = np.ones(np.shape(element_id)) # 1= active, 0 = inactive\n",
    "density = 2650 * np.ones(np.size(element_id))  # (kg/m3)\n",
    "abrasion_rate = 0 * np.ones(np.size(element_id)) # (mass loss /m)\n",
    "\n",
    "# Lognormal GSD\n",
    "medianD = 0.15 # m\n",
    "mu = np.log(medianD)\n",
    "sigma = np.log(2) #assume that D84 = sigma*D50\n",
    "np.random.seed(0)\n",
    "D = np.random.lognormal(\n",
    "    mu,\n",
    "    sigma,\n",
    "    np.shape(element_id)\n",
    ")  # (m) the diameter of grains in each parcel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to track sediment motion, we classify parcels as either active (representing mobile surface sediment) or inactive (immobile subsurface) during each timestep. The active parcels are the most recent parcels to arrive in the link. During a timestep, active parcels are transported downstream (increasing their `location_in_link`, which is a normalized value ranging from 0 to 1) according to a sediment transport formula. \n",
    "\n",
    "We begin by assigning each parcel an arbitrary (and small) arrival time and location in the link. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_arrival_in_link = np.random.rand(np.size(element_id), 1) \n",
    "location_in_link = np.random.rand(np.size(element_id), 1) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to the required parcel attributes listed above, you can designate optional parcel characteristics, depending on your needs. For example: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lithology = [\"quartzite\"] * np.size(element_id)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now collect the arrays into a dictionary of variables, some of which will be tracked through time (`[\"item_id\", \"time\"]`), and others of which will remain constant through time :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "variables = {\n",
    "    \"abrasion_rate\": ([\"item_id\"], abrasion_rate),\n",
    "    \"density\": ([\"item_id\"], density),\n",
    "    \"lithology\": ([\"item_id\"], lithology),\n",
    "    \"time_arrival_in_link\": ([\"item_id\", \"time\"], time_arrival_in_link),\n",
    "    \"active_layer\": ([\"item_id\", \"time\"], active_layer),\n",
    "    \"location_in_link\": ([\"item_id\", \"time\"], location_in_link),\n",
    "    \"D\": ([\"item_id\", \"time\"], D),\n",
    "    \"volume\": ([\"item_id\", \"time\"], volume)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With all of the required attributes collected, we can create the parcels DataRecord. Often, parcels will eventually transport off of the downstream-most link. To track these parcels, we have designated a \"`dummy_element`\" here, which has index value `-2`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "items = {\"grid_element\": \"link\", \"element_id\": element_id}\n",
    "\n",
    "parcels = DataRecord(\n",
    "    grid,\n",
    "    items=items,\n",
    "    time=[0.0],\n",
    "    data_vars=variables,\n",
    "    dummy_elements={\"link\": [NetworkSedimentTransporter.OUT_OF_NETWORK]},\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Run the NetworkSedimentTransporter\n",
    "\n",
    "With the parcels and grid set up, we can move on to setting up the model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timesteps = 20 # total number of timesteps\n",
    "dt = 60 * 60 * 24 *10 # length of timestep (seconds) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before running the NST, we need to determine flow direction on the grid (upstream and downstream for each link). To do so, we initalize and run a Landlab flow director component: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fd = FlowDirectorSteepest(grid, \"topographic__elevation\")\n",
    "fd.run_one_step()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we initialize the network sediment transporter: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nst = NetworkSedimentTransporter(    \n",
    "    grid,\n",
    "    parcels,\n",
    "    fd,\n",
    "    bed_porosity=0.3,\n",
    "    g=9.81,\n",
    "    fluid_density=1000,\n",
    "    transport_method=\"WilcockCrowe\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to run the model forward in time: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for t in range(0, (timesteps * dt), dt):\n",
    "    nst.run_one_step(dt)\n",
    "    \n",
    "    print(\"Model time: \", t/(60*60*24), \"days passed\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Plot the model results\n",
    "\n",
    "There are landlab plotting tools specific to the `NetworkSedimentTransporter`. In particular, `plot_network_and_parcels` creates a plan-view map of the network and parcels (represented as dots along the network). We can color both the parcels and the links by attributes.  \n",
    "\n",
    "Here, we demonstrate one example use of `plot_network_and_parcels`. For a thorough tutorial on the plotting tools, see [this notebook](../network_sediment_transporter/network_plotting_examples.ipynb).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can color links by values that we calculate. For example, if we are curious about the fate of sediment that started out on link 300 (in the far northwest corner of the basin), we might want to plot the total volume of sediment that originated on link 300 during a later timestep:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "timestep_of_interest = 6\n",
    "originating_link = 300\n",
    "\n",
    "#filter the parcels to calculate total volumes of only the parcels that originated in the chosen link\n",
    "parcelfilter = np.zeros_like(\n",
    "    parcels.dataset.element_id, dtype=bool\n",
    ")\n",
    "parcelfilter[:, timestep_of_interest] = (parcels.dataset.element_id[:,0] == originating_link)\n",
    "\n",
    "vol_orig_link = parcels.calc_aggregate_value(\n",
    "            np.sum,\n",
    "            \"volume\",\n",
    "            at=\"link\",\n",
    "            filter_array=parcelfilter,\n",
    "            fill_value=0.0\n",
    "        )\n",
    "\n",
    "fig = plot_network_and_parcels(\n",
    "    grid, parcels,\n",
    "    link_attribute=vol_orig_link, \n",
    "    link_attribute_title = \"Vol of sed originating on link x\",\n",
    "    network_linewidth = 5,\n",
    "    parcel_alpha = 0\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Non-network plotting\n",
    "The results of the NST can be visualized by directly accessing information about the grid, the parcels, and by accessing variables stored after the run of NST. \n",
    "\n",
    "As a simple example, we can plot the total volume of parcels on the grid through time. As parcels exit the grid, the total volume decreases.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "parcel_vol_on_grid = parcels.dataset[\"volume\"].values\n",
    "parcel_vol_on_grid[parcels.dataset[\"element_id\"].values==-2]=0\n",
    "\n",
    "#plt.figure(figsize=(8,6))\n",
    "plt.plot(np.asarray(parcels.time_coordinates)/(60*60*24), \n",
    "         np.sum(parcel_vol_on_grid, axis=0),\n",
    "         '-',\n",
    "         linewidth=3, \n",
    "         alpha=0.5\n",
    "        )\n",
    "\n",
    "plt.ylabel('Total volume of parcels on grid $[m^3]$')\n",
    "plt.xlabel('Time [days]')\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also plot individual parcel characteristics. The plot below shows the total transport distance of each parcel through the whole model run as a function of the parcel's grain size (during the final timestep). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.loglog(parcels.dataset.D[:,-1],\n",
    "         nst._distance_traveled_cumulative,\n",
    "         '.'\n",
    "        )\n",
    "plt.xlabel('Parcel grain size (m)')\n",
    "plt.ylabel('Cumulative parcel travel distance (m)')\n",
    "\n",
    "# Note: some of the smallest grain travel distances can exceed the length of the \n",
    "# grid by \"overshooting\" during a single timestep of high transport rate\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot below is an example of accessing variables associated with the grid (`grid.at_link.X`, or `grid.at_node.X`), as well as a variable associated with this instance of NetworkModelGrid (`nmg.X`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(grid.at_link[\"channel_slope\"],\n",
    "         nst.d_mean_active, \n",
    "         '.')\n",
    "plt.xlabel('Channel slope (m/m)')\n",
    "plt.ylabel('Mean grain size of active layer (m)')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
